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The phonon modes of the Frenkel-Kontorova model are 
studied both at the pinning transition as well as in the pinned 
(cantorus) phase. We focus on the minimal frequency of the 
phonon spectrum and the corresponding generalized eigen- 
function. Using an exact decimation scheme, the eigenfunc- 
tions are shown to have nontrivial scaling properties not only 
at the pinning transition point but also in the cantorus regime. 
Therefore the phonons defy localization and remain critical 
even where the associated area-preserving map has a posi- 
tive Lyapunov exponent. In this region, the critical scaling 
properties vary continuously and are described by a line of 
renormalization limit cycles. Interesting renormalization bi- 
furcation diagrams are obtained by monitoring the cycles as 
the parameters of the system are varied from an integrablc 
case to the anti-integrable limit. Both of these limits are de- 
scribed by a trivial decimation fixed point. Very surprisingly 
we find additional special parameter values in the cantorus 
regime where the renormalization limit cycle degenerates into 
the above trivial fixed point. At these "degeneracy points" the 
phonon hull is represented by an infinite series of step func- 
tions. This novel behavior persists in the extended version of 
the model containing two harmonics. Additional richnesses of 
this extended model are the one to two-hole transition line, 
characterized by a divergence in the renormalization cycles, 
nonexponentially localized phonons, and the preservation of 
critical behavior all the way upto the anti-integrable limit. 



I. INTRODUCTION 

The Frenkel-Kontorova (FK) model is one of the sim- 
plest prototypes for theoretical studies on incommensu- 
rate structures in solid-state physics |4|,[7],[l3| . The model 
consists of an infinite one-dimensional chain of balls in an 
external sinusoidal potential with an elastic interaction 
between nearest neighbors. If the mean distance between 
consecutive balls is not in rational proportion to the pe- 
riodicity of the sinusoidal potential, the corresponding 
state is called incommensurate. Incommensurate struc- 
tures appear in many physical systems such as quasicrys- 
tals, two-dimensional electron systems, magnetic super- 
lattices, and charge-density waves. 



Because of the well-known connection to the area- 
preserving standard map, the FK model is important also 
from the nonlinear-dynamics point of view. An incom- 
mensurate ground state of the FK model corresponds ei- 
ther to an invariant circle or a cantorus (invariant Cantor 
set) for the standard map. In the Kolmogorov-Arnol'd- 
Moser (KAM) phase, a configuration associated with a 
smooth circle can be described by an analytic hull func- 
tion which maps a rigid rotation orbit to the configu- 
ration. On the other hand, a cantorus corresponds to 
a purely discrete hull function representable as a con- 
vergent series of step functions jfj). Physically, a state 
with an analytic hull function is free to slide whereas a 
discrete hull function implies a non-zero minimum force 
(or field) required for depinning. The pinning transition 
corresponding to the breakup of an invariant circle is usu- 
ally refered as the transition by breaking of analyticity 
(TBA). 

The linear excitation or the phonon modes of the FK 
model describing the physical stability of the states are 
affected in an important way by the TBA. In the KAM 
phase and at the TBA, it is possible to find a transla- 
tional mode with zero frequency (phason mode) support- 
ing collective transport. The corresponding generalized 
phonon cigenfunctions arc "extended" (quasiperiodic) in 
the KAM phase and "critical" (neither extended nor lo- 
calized) with non-trivial scaling properties at the TBA. 
However, not much has been known about the eigen- 
functions in the cantorus regime where the frequencies of 
phonons are bounded away from zero. In a related quan- 
tum problem, describing an electron in a quasiperiodic 
potential, a similar transition results in localized eigen- 
functions [^0| . On the other hand, the absence of localiza- 
tion has been proven for a quasiperiodic potential taking 
a finite number of discrete values jlj]. Although the po- 
tential seen by the phonons is purely discrete in the can- 
torus regime, the potential still takes an infinite number 
of distinct values which leaves the question on the nature 
of the eigenfunctions open. There are recent speculations 
by Burkov et al. [Q on the possibility of having critical 
phonon modes at band edges of the phonon spectrum. 
However, in their numerical study, where the localiza- 
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tion character is infered by computing the participation 
ratio, it is rather difficult to make definite predictions. 

A great deal of effort has been put on studying the non- 
trivial scale invariance observed at the TB A |f^Jl^,|f^,^] . 
MacKay j2^] has adjusted his previous renormalization 
theory for the critical invariant circles to explain many 
scaling exponents in the neighborhood of the transition. 
However, it has not been clear how to formulate the 
scaling theory in the supercritical cantorus region where 
many relevant quantities like the actions or phase space 
scaling constants diverge in renormalizing plj . Sneddon 
et al. (2(| took a step in this direction by proposing a 
renormalization theory of the pinning threshold in the 
case of infinite viscosity. Working in the Fourier basis, 
they found three-term harmonic expansions for both the 
TBA renormalization fixed point and the new dynamic- 
threshold fixed point characterizing the depinning tran- 
sition at a critical field in the supercritical regime. The 
latter fixed point corresponded to vanishing elastic con- 
stants. 

In this paper, an exact decimation approach |l7| , ^9| , po| 
is used to study phonons in the FK model both at the 
TBA and beyond the transition. We show that the mini- 
mal frequency phonon mode remains critical throughout 
the cantorus regime. Furthermore, the scaling proper- 
ties of the corresponding generalized eigenfunctions vary 
continuously and are characterized by a line of limit cy- 
cles of the renormalization equations. We numerically 
obtain a full bifurcation diagram of these cycles as the 
system parameter is varied upto the anti-integrable limit 
where the elastic interaction disappears. The limit cycle 
degenerates to a trivial fixed point at the anti-integrable 
limit and at an infinite sequence of other parameter val- 
ues approaching the TBA from above. At this sequence 
of points, the phonon eigenfunctions are represented by 
infinite series of step functions and hence provide novel 
behavior in the study of systems with competing period- 
icities. The degeneracy points persist in the FK model 
with two harmonics and hence are believed to be generic 
in one-parameter families of similar models. 

We also investigate the scaling properties of one to two- 
hole transition ||10(] for cantori in a two-harmonic model. 
Using a multi-precision software package ^M , we show 
that the transition is signalled by a divergence in one 
member of the limit cycle of the renormalization tra- 
jectory. However, the finite members of the limit cycle 
vary continuously along the transition line in the two- 
parameter space. This transition line is distinguished 
from the rest of the cantorus phase by the fact that the 
phonon eigenfunctions vanish asymptotically. However, 
the localization is not exponential as suggested by van- 
ishing of the Lyapunov exponent. Another new feature 
existing only for a two-hole cantorus is the non-trivial 
anti-integrable limit for phonons, i.e. the phonons re- 
main critical all the way upto the anti-integrable limit. 

The FK model, its relation to the standard map and 



the phonon modes describing the physical stability of the 
FK model are briefly reviewed in Section II. Section III 
explains the decimation method which has been previ- 
ously applied to investigate the critical invariant circles 
in dissipative systems Jl7| and the self-similar wave func- 
tions in tight binding models [ ^9pC| ]. In Section IV we 
give an account of the numerical methods used in our 
study. The results of applying the decimation method 
to the FK model are reported in Section V. Analogous 
results on the extended FK model with two harmonics 
are found in Section VI. The main emphasis there is on 
the one to two-hole transition. The results and some 
open problems are discussed in Section VII. An appendix 
contains a detailed discussion on what happens to the 
phonons approaching the anti-integrable limit. 



II. GENERALITIES 

This section is about the FK model, terminology and 
basic results of Aubry's [HQ] theory. 

A configuration of the FK model is defined as the se- 
quence Xi, i € Z, where Xi gives the position of the ith 
ball. The energy of the configuration is given by 



where v is the generating function 
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v(xi,x i+ i) = -(x i+ i - Xi - a) +V(xi), 
k 

V(x) = — ^ cos(27ra;). 



(1) 

(2) 
(3) 



Because the external potential part V of the energy has 
the period unity, two configurations x and x can be iden- 
tified if there exist m, n £ Z such that the transformation 
{T mn x)i — Xi+m — n maps x on x. An equilibrium state 
with dW/dxi — for all i € Z is obtained by iterating 
the area-preserving standard map 

(x i+1 ,y i+1 ) = F(x l ,y l ) = (x. t + y i+ i,yi + V'(xi)), (4) 

where yi = Xi - x l -\. 

A minimum energy (m.e.) state is defined as a configu- 
ration x such that the energy for all other configurations 
which agree with x outside a finite range of i is greater 
than or equal to the energy of x. Every m.e. state is ro- 
tationally ordered which implies the existence of the 
mean spacing of the balls given by the rotation number 



(j = lim 



N 



(5) 



Incommensurate states correspond to irrational values of 
a. 

One is usually interested in recurrent m.e. states which 
are called ground states by Aubry There exists a 
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closed non-empty set of recurrent m.e. states of every 
mean spacing a £ R. If a is irrational, this set corre- 
sponds either to an invariant circle or a cantorus around 
the cylinder obtained by identifying x and x + 1. In 
other words, the FK model can have an infinite number 
of ground states all corresponding to the same invari- 
ant circle or cantorus of the standard map. An invariant 
circle can be described by a unique strictly increasing 
continuous hull function g : R — ► R so that 



Xi=g(i* + <f>), g(0 + l)=g(9) + h 



(6) 



where by varying the phase <f> one can choose the ground 
state one wants. There is a similar representation for a 
cantorus in terms of either the left- or right-continuous 
extension of a unique non-decreasing discontinuous hull 
function. 

The physical ground states are obtained by minimiz- 
ing the energy per ball over a. For the FK model, the 
parameter a can be tuned for each k so that the physical 
ground states have a desired mean spacing. 

The equation describing the physical stability of the 
ball configurations in the FK model is obtained by lin- 
earizing the dynamics d 2 Xi/dt 2 = —dW/dxi about an 
equilibrium configuration. The resulting equation for the 
small displacements £j is 

d% 
dt 2 



d 2 W 



dxidxj ^ 



(J) 



which can be written in the compact way d 2 £,/dt = 
— D 2 W£. The phonon modes are solutions of the form 
£(t) = exp(iujt)ip , where the (generalized) eigenvector tp 
of the operator D 2 W satisfies the equation 

ip i+1 + ^i-i - [2 + V"{xi)]i>i = -uj 2 ^. (8) 

The spectrum of D 2 W, given in terms of the squares 
of phonon frequencies to 2 , relates in an important way 
to the properties of the ground states of the FK model. 
First, it can be proven that the spectrum of D 2 W does 
not depend on the phase of an incommensurate ground 
state in the FK model f|. For an invariant circle, the 
infimum E m i n of the spectrum is zero whereas for a 
uniformly hyperbolic cantorus E min > H. In gen- 
eral, the configurations with F m - m > are physically at 
least metastable. Morover, uniform hyperbolicity implies 
structural stability, i.e. there cannot be any bifurcations 
of the cantorus at the corresponding parameter values. 
However, the cantorus can lose uniform hyperbolicity at 
a bifurcation where E m i n = 0. Such bifurcations appear 
in extended versions of the FK model containing two or 
more harmonics P, |l(i| , |l8| , p4| . 

III. PHONON DECIMATION 

The phonon equation (JsJ) can be interpreted as a dis- 
crete Schrodinger equation for the "wave function" ipi in 



the quasiperiodic potential 2 + V"(g(ia + (/)))■ Similar 
quasiperiodic linear-difference equations appear in many 
other nonlinear dynamics and solid state problems. We 
have developed a decimation method to deal with such 
equations |l^,|l^,^0| . The main advantage of the method 
is that the scaling properties of solutions can be described 
by bounded decimation functions. The decimation strat- 
egy is determined by the underlying "frequency" a. Our 
decimation can be implemented for any irrational a by 
choosing an arbitrary path in the Farey tree []l9|. How- 
ever, in this paper we will restrict ourselves to one of 
the noble rotation numbers, namely, a — 7 -2 , where 
7 = |(1 + V5) is the golden ratio. 

In the decimation scheme, all other lattice sites except 
those labelled by the Fibonacci numbers F n , F n+1 = F n + 
F n ^i (Fq = 0, F\ = 1), are decimated out beginning from 
an arbitrary initial site i. At the n th step, the decimated 
phonon equation has the form 



f n (i)ip(i + K+i) = f/>(i + F n ) + e n (i)i/)(i). 



(9) 



The phonon frequency is simply a parameter in this dec- 
imation equation which therefore should not be consid- 
ered as an eigenvalue problem. Using the additive prop- 
erty of the Fibonacci numbers, straightforward manupu- 
lations give the following exact recursion relations for the 
decimation functions e n and /„: 



fn+i(i) 



Ae n (i) 
l+Af n (i) 
f n -i{i + F n )f n (i + F n ) 



(10) 
(11) 



l + A/„(») 

A = e n _i(« + F n ) + f n -i(i + F n )e n (i + F n ). 

The above recursion relations can be iterated numerically 
for an arbitrary number of decimation steps provided the 
phonon frequency to and the sequence X4, which set the 
initial functions e.i(i) and fii^)-, can be calculated with 
arbitrary precision. 

There are a number of reasons why the iteration of the 
recursion relations ( p~0| - |TT| ) is superior to the direct itera- 
tion of the phonon equation (^) : 

1) Often the iteration of the recursion relations is numer- 
ically more stable than that of Eq. (||). 

2) The possible self-similarity observed by monitoring the 
behaviour of ipi over the range i £ [— F n , F n ] with increas- 
ing n can be captured by a simple asymptotic limit cycle 
for the decimation functions e n and /„. 

3) An unknown or inaccurate parameter in the equation 
(^) and the asymptotic limit cycle can be simultaneouly 
determined self-consistently (see the next section for fur- 
ther details). In other words, the decimation equations 
themselves provide a new method to compute e.g. eigen- 
values or phase boundaries upto machine precision. 

The decimation scheme can also be use in describ- 
ing the self-similar fluctuations of exponentially localized 
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modes. Here the localized eigenfunction is written as 



--y\i\ 



r\i which results in the following equation for 



"Trji+i + e^-i - [2 + V"(xi)]TH = -u\ 



(12) 



The above equation can be studied using the decimation 
scheme with the inverse localization length 7 being an 
additional parameter. Again 7 can be determined self- 
consistently along with the limit cycle. Q 

Although the above scheme with a discrete lattice in- 
dex i is sufficient for the present purpose, we would like 
to point out that an analogous form of the recursion re- 
lations in terms of a continuous variable is needed if one 
wants to solve renormalization limit cycles by the Newton 
method MM. 



IV. NUMERICAL METHODS 

In this paper, the properties of the solutions to Eq. 
(^) are studied at the minimal frequency u> — \[E min . 
In the KAM phase, the ground state configurations of 
the model exhibit a zero frequency phonon. In this case, 
the corresponding eigenfunction can be determined very 
accurately. At the TBA, the precision in obtaining the 
asymptotic scaling properties is limited by the precision 
of the critical k value and the sequence Xi . On the other 
hand, the iteration of the decimation equations for an 
arbitrary value of the nonlinearity parameter k in the 
cantorus regime is mainly limited by the precision of the 
phonon frequency. However, as pointed out in the previ- 
ous section, the existence of asymptotic limit cycles for 
the decimation functions can be used to obtain the fre- 
quency to high precision so that the remaining source of 
numerical error is the inaccuracy of the ball positions. 
This error can be made very small by approximating a 
ground state with very long orbit. 

There are essentially two steps involved in obtaining 
the scaling properties of the phonon eigenmodes. The 
first step is to obtain an estimate for the minimal fre- 
quency of the phonon spectrum using an orbit of moder- 
ate size (e.g. 610) in the diagonalization procedure. The 
frequency estimate has to be precise enough in order to 
allow so many decimation steps that the asymptotic cy- 
cle is already approached. Usually the diagonalization 
procedure is the most CPU time consuming part in the 
numerical calculations. The second part is to approxi- 
mate the cantorus with a very long orbit (e.g. of the 
length 46368) and to improve the frequency estimate to 
very high precision so that the recursion relations ( fLOf - 
pT| ) can be iterated many times to obtain the asymptotic 
cycle with several significant digits. 

We follow Baesens and MacKay (see Appendix A of ref. 
pfj|) in calculating an initial estimate for E m i n and the 
sequence Xi for the decimation equation (H) . An invariant 



circle or a cantorus of rotation number 7 -2 is approxi- 
mated by periodic orbits of rotation numbers F n -2/F n . 
It is sufficient to consider symmetric periodic orbits which 
fall into four classes |U^: 1) even period with xq = 0, 2) 
even period with xq = —x\, 3) odd period with xo = 0, 
and 4) odd period with x = 1/2. These classes are use- 
ful in obtaining initial conditions for the iteration of Eq. 
(||) as the phonon modes in each class are either sym- 
metric or antisymmetric [|l0|. In the KAM phase and at 
the TBA, the invariant circle describing the ball config- 
urations can be approximated by any of the above four 
classes. However, a cantorus can be approximated only 
by class-2 (or class-4 if the cantorus has only one hole) 
orbits. The classes 1 and 3 (and 4 in the region of two 
holes, see section VI) give some homoclinics in addition to 
the cantorus so E m i n can be different for them. Baesens 
and MacKay explain how to obtain the orbit Lyapunov 
exponent or the residue needed in locating the TBA and 
the one to two-hole transition on the parameter axis. It 
should be noted that in the cantorus phase, the standard 
map has a positive Lyapunov exponent which near anti- 
integrable limit exhibits a logarithmic divergence with 
the nonlinearity parameter k. 

One of our main tasks was to decide on whether the 
eigenfunctions were exponentially localized or not. Thou- 
less has derived a formula for the inverse localization 
length [^5). However, it practice the estimates obtained 
with this formula remain rather crude as it requires con- 
sidering differences between all eigenvalues. Initially, we 
carried out the decimation varying the value of the in- 
verse localization length 7 but we found that only 7 = 
showed convergence towards a non-trivial asymptotic p- 
cycle for the decimation functions at lu = \f~E min . As- 
suming 7 = 0, the value of E m i n was then sequentially 
improved with the secant method by requiring that the 
decimation functions asymptotically converge to a period 
p limit cycle. In particular, we used / n (0) = f n - P {0) with 
e.g. n — 17, 20, 23, where the cantorus was approximated 
with an orbit of the length F23 =46368. Table I compares 
E m i n obtained in this way at k = 2 with diagonalization 
results. One should note that the above self-consistent 
calculation of E min is not only more accurate but also 
takes only a small fraction of the CPU time needed for 
the diagonalization of a matrix of the size 10946/2. All 
the computations were done using quadruple precision. 

The critical line for the one to two-hole transition was 
determined by the symmetry breaking bifurcations as in 
ref. 1 1 1 . Table II lists the bifurcation points of periodic 
orbits along the line k^ — 4fci . It can be seen in this table 
that due to superexponential convergence, all decimals of 
the quadruple precision are already exhausted at the or- 
bit length 377. This size is highly inadequate because it 
limits the maximum number of decimation iterations to 
be 14. With 14 iterations, we are still in the transient 
regime and the existence of a limit cycle for the decima- 
tion functions is far from obvious. Moreover, the nature 
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of phason eigenfunction or the fate of the orbit Lyapunov 
exponent cannot be decided on such a short orbit. To 
overcome this problem, we used a multi-precision soft- 
ware package JTl| to compute the orbit and the phonon 
eigenfunction along the one to two-hole transition line 
(see Section VI for further details). 



V. RENORMALIZATION BIFURCATION 
DIAGRAM 

The two-step numerical procedure discussed in the pre- 
vious section was applied to study the scaling properties 
of the phonon eigenfunctions for a multitude of values of 
k starting from the KAM phase upto the anti-integrable 
limit k — * oo. This limit corresponds to vanishing elas- 
tic interaction H which implies that for a m.e. state 
all balls lie at the potential energy minimum x = 1/2. 
By introducing the parameter r — 2 arctan(fc)/7r one can 
map the whole fc-axis to the r-interval [0, 1], where r = 
and r = 1 are the integrable and anti-integrable limits, 
respectively. 

In the KAM region, where the invariant circle is 
smooth, one obtains the same trivial fixed point for 
all classes of approximating orbits. It can be easily 
solved from the fixed point equation: f n (i) = & and 
e n (i) = —a 2 . It describes an extended quasiperiodic 
phonon eigenfunction. 

In general, choosing the class fixes the phase (j> in Eq. 
(|J) and thus different classes can lead to different deci- 
mation cycles. This is the case at the TBA, where the 
class 1 and 3 orbits having xq on the dominant symmetry 
line lead to a nontrivial universal fixed point M] whereas 
the class 2 and 4 orbits correspond to 3-cycles. In Fig. 
1 we plot the critical phason (u> = 0) mode ipi obtained 
by approximating the invariant circle by the class 1 or 3 
orbits which causes the main peak to be located at i = 0. 
The non-trivial fixed point and the corresponding self- 
similarity clearly distinguishes the TBA from the KAM 
regime. 

In the cantorus regime, both classes 2 and 4 lead to 
the same decimation 6-cycle for fixed k. Fig. 2(a) shows 
the variation of the 6-cycle as we change the parameter r 
from the KAM region to the anti-integrable limit. An in- 
teresting feature of the diagram is the infinite number of 
loops separated by degeneracy points where the decima- 
tion functions approach asymptotically the same trivial 
limit cycle as in the KAM phase. Coming down from the 
anti-integrable limit, the first such degeneracy is found 
at r = 0.67052674.... Although the decimation func- 
tions come close to the trivial fixed point already near 
r = 0.84, they behave parabolically in this region forming 
an avoided crossing. As shown in the blowup of Fig. 2(b), 
the behaviour is very different around r = 0.67052674 
where the curves behave linearly as a function of the pa- 
rameter r all passing through the inverse golden mean 



at the same point. Table III gives a list of degeneracy 
points. We have checked that they all correspond to a 
real crossing of the curves in the bifurcation diagram. 
We conjecture that there are an infinite number of such 
points accumulating at the critical parameter value for 
the TBA. 

As can be seen in Table III, the separation between 
subsequent degeneracy points appear to asymptotically 
scale by the universal scaling constant 6 = 1.62795... for 
the TBA 21,22]]. This suggests that the observed loop 
pattern is related to the renormalization dynamics along 
the unstable manifold of the TBA fixed point. This in 
turn causes one to think that there should be some prop- 
erty of the configurations which makes the degeneracy 
points special. If such a property exists, it has to be 
rather subtle as the orbit hull function is given by quali- 
tatively the same kind of series of step functions through 
the whole cantorus region. 

An important property associated with the degener- 
acy points is revealed by plotting the phonon hull, i.e. ipi 
versus {ia} where {icr} is the fractional part of ia. Fig. 
3(a) shows this plot at one of the degeneracy points. The 
resulting graph appears to be represented by an infinite 
series of step functions. Similar plots are obtained at 
other degeneracy points. In other words, the phonon hull 
shows similar characteristics as the orbit hull function g 
at these points. It is interesting to note that this is not 
only true at the degeneracy points but also in the KAM 
phase where both the orbit and phonon hull functions 
are smooth. This suggests that trivial asymptotic dec- 
imation functions in general signify qualitatively equiv- 
alent behavior in the quasiperiodic "potential" term of 
the linear difference equation and in the corresponding 
(generalized) eigenfunction. 

In the cantorus regime outside the degeneracy points, 
the phonon hull consists of a fractal set of points resem- 
bling the phonon hull function at the TBA. Fig. 3(b) 
shows how the character of the phonon hull changes as 
we move away from the degeneracy points. The exis- 
tence of a decimation limit cycle implies self-similarity of 
the phonon eigenfunction. This self-similarity is clearly 
seen after adjusting the phase factor so that the eigen- 
functions remain bounded asymptotically ||. Interest- 
ingly, the resulting decimation trajectories converge on 
a 4-cycle. The 4-cycle appears to be more fundamental 
(and interesting because it breaks the number theoretic 
even-odd 3-cycle) than the 6-cycle. It means there is a 
point on the cantorus around which the phonon mode 
is self-similar although none of the approximating orbits 
pass through it. However, since the qualitative features 
of the bifurcation diagram (loops and degeneracy points) 
are similar using the 4 or 6-cycle, for simplicity, our bi- 
furcation diagrams are shown with 6-cycle. 
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VI. TWO-HARMONIC MODEL AND THE ONE 
TO TWO-HOLE TRANSITION 

Next we study the phonons and the associated renor- 
malization bifurcation diagram in the two-harmonic FK 
model with potential 

V( x ) — cos(27ra;) + cos(47ra). (13) 

Our motivation for this study is two-fold: firstly, to see if 
loops and degeneracy points are generic in such models 
and secondly, to obtain the renormalization bifurcation 
diagram near the one to two-hole transition for cantori 
studied recently by Baesens and MacKay fllpf . 

As stated in their paper, a gap in a cantorus is a pair of 
distinct points of the cantorus whose forward and back- 
ward orbits converge together. An orbit of gaps is called a 
hole. In the standard map with single well potential, the 
cantorus is uniquely determined by the rotation number 
(if there is no invariant circle with that rotation num- 
ber) and contains only one hole. However, in the two- 
harmonic model the cantori for a given rotation num- 
ber are parametrized by a £ [0, 1] near a non-degenerate 
double- well anti-integrable limit with | fci| < \k 2 \ ||. Here 
a gives the fraction of the natural measure of the cantorus 
in one of the wells. In other words, there are uncount- 
ably many cantori of each rotation number. They have 
generically two holes although there is a dense set of val- 
ues of a corresponding to only one hole. Using symmetry 
properties it can be argued for k 2 > that the m.e. can- 
torus has an equal fraction of points at both wells near 
the anti-integrable limit tj. Thus Baesens and MacKay 
fiof focus on a = 1/2. Numerically they observe a criti- 
cal line in the k\ — k 2 parameter space corresponding to 
the one to two-hole transition for the symmetric cantorus 
with E min = 0. 

Fig. 4(a) shows the renormalization bifurcation dia- 
gram obtained by varying the parameter 

r = — arctan \jk\ + k\ (14) 

along the line k 2 = 4fci . We followed the same numerical 
methods as in the previous section and approximated the 
cantorus by class-2 orbits (the class-4 orbits give a ho- 
moclinic in addition to the cantorus in the region of two 
holes). The figure clearly shows the existence of loops 
and degeneracy points near the TBA thus establishing 
the result that their existence is generic in one-parameter 
systems. 

The transition from one to two holes is signalled by a 
divergence in one of the six members of the limit cycle. 
The other members of the limit cycle are found to re- 
main finite. Furthermore, the bifurcation diagram in the 
neighborhood of the one to two-hole transition exhibit a 
rather strange pattern due to the fact that various mem- 
bers of the limit cycle cross each other as the parameter r 



is varied. However, the limit cycle retains the period six 
at these crossing points. The crossing pattern becomes 
more and more complicated near the 1-2 hole transition. 
There can be a multitude of crossings in a very narrow pa- 
rameter interval which makes the behaviour of the curves 
look discontinuous at some points. However, blowups re- 
veal that the curves are in fact smooth (see Fig. 4(b)). 

We have also studied the critical line corresponding to 
the hole transition in the k\ — k 2 space. Our numerical 
results strengthen the previous conjecture |l(J that the 
Lyapunov exponent of the extended standard map van- 
ishes along the hole transition (see Table IV). The deci- 
mation functions obtained with class 2 orbits converge on 
a 6-cycle with one member of the cycle diverging while 
the remaing elements of the 6-cycle remain bounded and 
vary continuously along the transition line. Numerically 
the phason mode is found to be peaked around the lat- 
tice site i c corresponding to Xi c which is closest to 1/2. 
Therefore, on the transition line it is natural to look at 
class-4 orbits because they have the 0th point exactly at 
1/2. However, we find no non-trivial limit cycle with the 
class-4 orbits reflecting the fact that the phason mode 
vanishes asymptotically as i — > ±oo. This localization 
is non-exponential with the inverse localization length 
7 = 0. In Fig. 5 we plot the phason mode in both 
linear-log and log-log scales. These plots show that the 
localization is slower than exponential but faster than 
geometric. To our knowledge, this peculiar behaviour of 
the eigenfunctions has not been observed before. 

Another unique aspect of the two-harmonic model is 
the fact that it is possible to show the existence of crit- 
ical phonon modes in the anti-integrable limit. This is 
clearly seen in Fig. 4(a) where the phonon eigenfunctions 
are characterized by a non-trivial limit cycle even ap- 
proaching r = 1 . As explained in Appendix, for a generic 
symmetric two- well potential the universality class in the 
anti-integrable limit is determined by a, a, and two addi- 
tional parameters setting the location of the minima and 
the ratio of the third and second derivative of the poten- 
tial at the minima. For the potential ( |l3| ) the latter two 
parameters are both fixed by the ratio k2/kl. 

VII. DISCUSSION 

Based on exact renormalization methods which can be 
implemented to very high precision on the computer, we 
have shown that the phonon eigenfunction at the mini- 
mal frequency is self-similar with scaling properties that 
vary continuously within the cantorus regime. A novel 
feature of our study is the existence of a sequence of 
points where the phonon eigenfunctions are represented 
by series of step functions. An important aspect of the 
work described here is the new method to compute eigen- 
values. The determination of the bifurcation diagram 
near the TBA was a challenging numerical problem and 



G 



would have been impossible without the self-consistent 
approach provided by the decimation scheme. 

It is interesting to compare our results with those for 
a related electron problem The TBA in the FK 

model is somewhat analogous to the metal-insulator tran- 
sition for two-dimensional electrons due to the formal 
similarity between the phonon equation and the corre- 
sponding equation for the electron wave functions. How- 
ever, the innocent looking difference in the quasiperi- 
odic potential leads to significantly different behaviour 
beyond the transition. In the "insulating" phase of the 
electron problem, the wave functions are exponentially 
decaying with self-similar fluctuations characterized by 
a unique renormalization fixed point. In contrast, the 
phonon eigenfunctions in the cantorus regime exhibit self- 
similarity without exponentially decaying envelope and 
can be described by a line of renormalization limit cy- 
cles. Furthermore, the cantorus phase contains an in- 
finite number of special points where the renormaliza- 
tion flow converges to the same trivial fixed point as in 
the KAM phase or in the "metallic" phase of the elec- 
tron problem. The anti-integrable limit is trivial for the 
standard FK model but non-trivial in the region of a 
two-hole cantorus. In comparison, the counterpart of the 
anti-integrable limit in the electron problem, the strong 
coupling limit is always highly nontrivial. 

The above described degeneracy points remain one of 
the main mysteries of our discoveries. In particular, why 
do the degeneracy points appear in the cantorus phase 
starting arbitrarily close to the TBA and why do they dis- 
appear beyond a certain value of the nonlinearity param- 
eter? We have tried to correlate them with some special 
characteristics of the ball configurations of the FK model 
but without success. Furthermore, the complexity of the 
renormalization flow near the one to two-hole transition 
is very puzzling and devoids of any explanation. 

Our results have significance in the theory of an almost 
periodic eigenvalue problem. Degeneracy points are the 
first example in this class of problems whose solutions 
can be represented by an infinite series of step functions. 
Another curious feature are the non-exponentially decay- 
ing eigenfunctions at the one to two-hole transition in an 
extended FK model. 

In this paper, we have explored only the parameter re- 
gion &2 > of the extended FK model. However, there is 
a really rich pattern of cantorus bifurcations for &2 < 
|],[l8|. In particular, Schellnhuber et al. have re- 
ported exponential localization of the phason mode at a 
boundary of metastability in an extended version of the 
FK model containing three harmonics. The metastabil- 
ity is related to vanishing of the gap of the phonon spec- 
trum which in that case signals annihilation of the can- 
torus configuration. It would be interesting to apply the 
decimation method at this and other similar bifurcations 
and determine the scaling properties of the exponentially 
localized eigenfunctions. 



One of the most intriguing features of studying the 
phonons in the FK model is that various different sub- 
jects like condensed matter physics, nonlinear dynamics, 
almost periodic Schrodinger operators, and renormaliza- 
tion theory are drawn together. We hope that the results 
of this paper motivate further research in each of these 
fields. 
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APPENDIX: PHONONS IN THE 
ANTI-INTEGRABLE LIMIT 

Here we explain the origin of the non-trivial critical be- 
havior of the phonons for the double-well anti-integrable 
limit (AIL). This will account for the differences at r = 1 
in Fig. 2(a) and the Fig. 4(a). In particular, we derive 
an equation for the critical phonons corresponding to a 
two- hole cantorus approaching the AIL. 

The key point of the derivation is the fact that the 
limiting form of the potential seen by the phonons ap- 
proaching the AIL is determined not by the locations of 
the balls at the AIL (which gives a constant contribution) 
but by the small deviation in the ball positions from those 
at the AIL. 

It is useful to multiply the phonon equation (||) by the 
factor 



t ■■ 



I 



y/kf + k 



(Al) 



so that all terms in the equation remain bounded in the 
AIL t -> 0: 



+ tipi-i - [2t + U"{xi)]tl)i = -u} 2 ipi 



(A2) 



where U(x) — tV(x) and a) 2 = tco 2 . It is important that 
U does not depend on the value of t if the ratio k^/kx 
is fixed. We now expand U(x) around the ball position 
at the AIL which we denote as e$. Since U'(ei) — 0, we 
have 



U'(x i ) = U"(e i )Z i + 0(g), 



(A3) 



where = (a;, — e,) is the small deviation in the ball 
configuration from that of the AIL. Using the equation 



U'(xi) = t(x i+ i - 2x l + Xi-i) 



(A4) 
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we notice that £j changes linearly with t near the AIL. 
Therefore, upto the first order in t, we have 



tdi 



(A5) 



where di — e 2 +i — 2e^ + ei-\. The next step is to put this 
together with the expansion 



U"(x i ) = U"(e i ) + U'"(e i )Z i +0(g) 



(A6) 



into Eq. (|Aj). Because U"(ei) is equal for all i if fc 2 > 0, 
it is possible to match the zeroth and first order terms in 
t so that one obtains the asymptotic equation 



i>i+i + + V A i L (i)i>i = (2 + C)ipi, 
where C is the limit of — oj 2 + U"(ei)/t as t - 

U"'(ei)di 



V AIL (i) 



U"(ei) 



(A7) 
and 

(A8) 



In the case of a single well, the potential is even around 
the minimum and the resulting asymptotic phonon equa- 
tion becomes simply 



VA:+i +ipi-i = (2 + C)Vi 



(A9) 



E m in corresponds to the supremum of the spectrum 
which in this case is 2, i.e. C = 0. 

In the case of two wells for \k\\ < k^ a straightforward 
calculation shows that 



±- 



Gtt 



(A10) 



where the + (— ) sign corresponds to {et\ > 1/2 ({ej} < 
1/2) . Let us next study the possible values of di assuming 
that a = (3 — \/5)/2. The hull function which gives the 
set a at the AIL reads |j|jic]] 



9e(0) = 



1 



Int(6),0< {9} < a 
m + lnt(0),a< {9} < 1 ' 



(All) 



where x m is the minimum in (0,1/2) satisfying 
cos(27ra; m ) = — kx/k^- In the symmetric case a = 1/2 
it is easy to see that di can have only the values 

±2x m , 0i G {(0, 1/2 - a), (1/2 + a, 1)} 
±{Ax m - l),6i G {(1/2 - a, a), (1 - a, 1/2 + a)}} 

(A12) 

where #j = ia + <p and the ± sign cancels the one com- 
ing from Eq. (AlC). Because the potential Vail takes 
only three different values, the generalized eigenfunctions 
cannot be localized p3 . E m i n corresponds to the supre- 



mum of the spectrum of Eq. (A7). The class 2 orbit 
corresponds to choosing the phase = — a/2 in the hull 



function flj0| . Decimation with this phase leads to the 
non-trivial 6-cycle found by extrapolating the curves of 
a bifurcation diagram to the AIL. 

It is interesting to note that Vail is determined solely 
by cr, a, x m , and a multiplying factor. Therefore, other 
similar two-well potentials for which these parameter 
agree would give rise to the same universality class in 
the AIL. 
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FIG. 3. Phonon hull at the second degeneracy point 
= 0.62221... (a) and away from it at r = 0.62 (b). 



FIG. 4. Renormalization bifurcation diagram for the ex- 
tended FK-model with ki = 4fci. (a) shows the elements of 
the limit cycle that remain bounded by 2 near the one to 
two-hole transition and (b) shows an example of how com- 
plicatedly the curves can cross one another in a very narrow 
parameter interval. For an analogous configuration bifurca- 
tion diagram see Fig. 1 of ref. 10. 



FIG. 5. Phason eigenfunction with xo = 1/2 at the one to 
two-hole transition plotted in linear-log (a) and log-log (b) 
scales. These plots show that the localization is between ex- 
ponential and power-law. 



TABLE I. Estimates for the infimum of the phonon spec- 
trum E m in at k = 2. The values in the third column were 
obtained by diagonalizing a matrix of size F n /2 as described 
in ref. 10. The values in the last column were determined 
by requiring that / n (0) = /n-e(0) where the sequence Xi was 
approximated by the class-2 orbit of the length 46368. 



(diag.) 



(dec.) 



14 
17 
20 
23 



610 

2584 
10946 
46368 



0.5935392137880 
0.5935394985221 
0.5935394812326 



0.5935394324321 
0.5935394834426 
0.5935394820297 
0.5935394820545 



FIG. 1. Universal phason mode at the TBA with xo = 0. 



FIG. 2. (a) Renormalization bifurcation diagram showing 
the asymptotic behaviour of the decimation function / as a 
function of the parameter r — 2 arctan(fc)/-7r. At each pa- 
rameter value, / n (0) is plotted for at least (p + 2) different 
decimation levels n which shows the accuracy of the limit 
cycle (of length p) obtained. The 3-cycle at the TBA is plot- 
ted with darker points. The curves behave parabolically near 
r = 0.84 and do not intersect there, (b) Blowup around the 
first degeneracy point where the branches of the 6-cycle meet 
at an intersection point. 
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TABLE II. Symmetry breaking bifurcations of class 4 or- 
bits along the line &2 = 4fci. The bifurcation points have 
been calculated using multiprecision software based on dou- 
bly precision arithmetics. We have checked consistence with 
a direct calculation using quadruple precision. The bifurca- 
tion points converge fast to the critical parameter value for 
the one to two-hole transition. Note that the corresponding 
value 0.429439737992501223898752836653609 obtained in ref. 
10 appears to be accurate only upto double precision. 



F„ ki 



377 0.4294397379925012239000251474523975392588953 
5245566429536, 

987 0.4294397379925012239000251474523975392588953 
7931282057125515563523777678 
714122115567680849012510667905, 

1597 0.4294397379925012239000251474523975392588953 
7931282057125515563523777678 
9131402934114670775341029280082904661719, 



TABLE III. Degeneracy points approaching the critical pa- 
rameter value r c from above and their scaling. 



n r n (r„_i - r c )/(r n - r c ) 

1.00000000 

1 0.67052674 2.8336 

2 0.62221221 1.3678 

3 0.56379391 1.8008 

4 0.53923590 1.5075 

5 0.51909955 1.7126 

6 0.50876597 1.5765 

7 0.50162025 1.6630 

8 0.49754645 1.6076 



TABLE IV. Lyapunov exponent A of the F n - 2 /F n orbit 
(class 2) at the one to two- hole transition point with = 4fci . 
The transition point was approximated by the last parameter 
value in Table II. The symmetry breaking points of the odd 
cycles converge faster than those of the even cycles which jus- 
tifies the inclusion of the Lyapunov exponent of the 2584-orbit 
in this table. 



F n A 

34 0.38572449 

144 0.20463797 

610 0.10618378 

2584 0.05450141 
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